%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 4 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: due to the use of a finite number of simulations 
%%% to evaluate (k1,k2,k3), the final results are random 
%%% and will not correspond exactly to those in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Define variables
nobs = 10^4;%We use nobs=10^7 to create Figure 4 in the paper
nu1 = 4;
nu2 = 8;
T = 100;
maxN = 60;
M = eye(T)-1/T;

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%% Initialize vectors
Expk1nu4 = zeros(maxN-1,1);
Expk2nu4 = zeros(maxN-1,1);
Expk3nu4 = zeros(maxN-1,1);
Expk1nu8 = zeros(maxN-1,1);
Expk2nu8 = zeros(maxN-1,1);
Expk3nu8 = zeros(maxN-1,1);
k1nu4 = zeros(maxN-1,1);
k2nu4 = zeros(maxN-1,1);
k3nu4 = zeros(maxN-1,1);
k1nu8 = zeros(maxN-1,1);
k2nu8 = zeros(maxN-1,1);
k3nu8 = zeros(maxN-1,1);

%%% Loop on nobs random draws to evaluate the expectations in (k1,k2,k3)
for n = 1:nobs
    if rem(n,1000)==0
       n
    end
    Y = randn(T,maxN);
    lam1 = sqrt((nu1-2)./chi2rnd(nu1,T,1));
    lam2 = sqrt((nu2-2)./chi2rnd(nu2,T,1));
    lam1Y = lam1.*Y; 
    lam2Y = lam2.*Y; 
    slam1Y = sum(lam1Y); 
    slam2Y = sum(lam2Y); 
    for j=2:maxN
        %\nu = 4
        lamYj = lam1Y(:,1:j);
        Wn = inv(lamYj'*M*lamYj);
        Wn2 = Wn^2;
        slamYj = slam1Y(1:j);     
        Expk1nu4(j-1,1) = Expk1nu4(j-1,1) + trace(Wn)/nobs;
        Expk2nu4(j-1,1) = Expk2nu4(j-1,1) + trace(Wn2)/nobs;
        Expk3nu4(j-1,1) = Expk3nu4(j-1,1) + slamYj*Wn2*slamYj'/nobs;
        %\nu = 8
        lamYj = lam2Y(:,1:j);
        Wn = inv(lamYj'*M*lamYj);
        Wn2 = Wn^2;
        slamYj = slam2Y(1:j);
        Expk1nu8(j-1,1) = Expk1nu8(j-1,1) + trace(Wn)/nobs;
        Expk2nu8(j-1,1) = Expk2nu8(j-1,1) + trace(Wn2)/nobs;
        Expk3nu8(j-1,1) = Expk3nu8(j-1,1) + slamYj*Wn2*slamYj'/nobs;
    end
end

%%% Compute (k1,k2,k3)
for j = 2:maxN
    N = j;
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    rho = N/T;
    %\nu = 4
    y = @(x) (nu1-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu1/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu1/(nu1-2)));
    varphi = 2*eta^2*(1-rho)/(nu1-eta*(nu1-2));
    k1nu4(j-1,1) = Expk1nu4(j-1,1)*a1/eta;
    k2nu4(j-1,1) = Expk2nu4(j-1,1)*a2/varphi;
    k3nu4(j-1,1) = Expk3nu4(j-1,1)*a3/eta;
    %\nu = 8
    y = @(x) (nu2-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu2/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu2/(nu2-2)));
    varphi = 2*eta^2*(1-rho)/(nu2-eta*(nu2-2));
    k1nu8(j-1,1) = Expk1nu8(j-1,1)*a1/eta;
    k2nu8(j-1,1) = Expk2nu8(j-1,1)*a2/varphi;
    k3nu8(j-1,1) = Expk3nu8(j-1,1)*a3/eta;
end

%%% Create Figure 4
Nvec = 2:maxN;
fig = figure;
%Subplot 1
subplot(2,2,1)
plot(Nvec,k1nu4,'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(Nvec,k2nu4,'--b','LineWidth',2);
plot(Nvec,k3nu4,':r','LineWidth',2);
hold off
legend('$k_1/\eta$','$k_2/\varphi$','$k_3/\varphi$','Location','northeast',...
       'interpreter','latex');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=4$','interpreter','latex');
xlabel('$N$','interpreter','latex');
xlim([2 60]);
ylim([0.995 1.08]);
xticks([2 10 20 30 40 50 60]);
yticks([1 1.02 1.04 1.06 1.08]);
%Subplot 2
subplot(2,2,2)
plot(Nvec,k1nu8,'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(Nvec,k2nu8,'--b','LineWidth',2);
plot(Nvec,k3nu8,':r','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\nu=8$','interpreter','latex');
xlabel('$N$','interpreter','latex');
xlim([2 60]);
ylim([0.9975 1.03]);
xticks([2 10 20 30 40 50 60]);
yticks([1 1.01 1.02 1.03]);
fontsize(fig,16,"points");